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Abstract 

Stochastic quantization offers the opportunity to simulate field theories 
with a complex action. In some theories unstable trajectories are preva- 
lent when a constant stepsize is employed. We construct algorithms for 
generating an adaptive stepsize in complex Langevin simulations and find 
that unstable trajectories are completely eliminated. To illustrate the gen- 
erality of the approach, we apply it to the three-dimensional XY model at 
nonzero chemical potential and the heavy dense limit of QCD. 
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1 Introduction 



Nonperturbative simulations of field theories with a complex action are difficult, 
since importance sampling based techniques break down. Stochastic quantisation 
and complex Langevin dynamics (TJ |2j |3l H] can potentially evade this problem, as 
it is not based on a probability interpretation of the weight in the path integral. 
Among others, this is relevant for theories with a sign problem due to a nonzero 
chemical potential El El [5] and for the dynamics of quantum fields out of 
equilibrium QUI ED]- Other recent activity can be found in Refs. [T21 fT3] . 

As is well known since the 80s, there are a number of problems associated 
with complex Langevin dynamics, see e.g. Refs. p3j [T5| [16]. These can roughly 
be divided under two headings: instabilities and convergence. The first problem 
concerns the appearance of instabilities when solving the discretized Langevin 
equations numerically. Sometimes, but not always, this can be controlled by 
choosing a small enough stepsize. The second problem pertains to convergence. 
In some cases complex Langevin simulations appear to converge but not to the 
correct answer (see e.g. Ref. [E]). In order to disentangle these issues, we tackle 
in this paper the first one and present adaptive stepsize algorithms that lead to a 
stable evolution and are not constrained to very small stepsizes only. A discussion 
of the second problem is deferred to future publications. 

The paper is organized as follows. In Sec. |2] we introduce the problem, in- 
dicate why an adaptive stepsize might be necessary and outline the basic idea 
behind the algorithms, extending the ideas of Refs. [T5j IT6] . To avoid notational 
cluttering we use a real scalar field, but we emphasize that the method is more 
generally applicable. We then present two algorithms implementing the basic 
idea, and apply them to the three-dimensional XY model at finite chemical po- 
tential in Sec. |3] and the heavy dense limit of QCD in Sec. HJ The latter has 
previously been studied with complex Langevin dynamics in Ref. [6]. We show a 
few selected results to indicate the applicability of the approach. As mentioned 
above, convergence will be discussed elsewhere. 

2 Adaptive stepsize 

Consider a real scalar field <fi with the Langevin equation of motion 



Here d is the supplementary Langevin time, —5S/ 5<j) is the drift term derived from 
the action S, and i] is Gaussian noise. The fundamental assertion of stochastic 
quantisation is that in the infinite (Langevin) time limit, noise averages of ob- 
servables become equal to quantum expectation values, defined via the standard 
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Figure 1: Example of a classical flow diagram in the XY model at nonzero chem- 
ical potential (/i = 2). The arrows denote the normalized drift terms (K R ,K l ) 
at ((f^,^ 1 ). The dots are classical fixed points. 



functional integral, 



lim (OWM,, - 1 Ziyi (22) 
■9-¥oo 1 J D(p e 



where the brackets on the left denote a noise average. 

If the action is complex the drift term becomes complex and so the field 
acquires an imaginary part (even if initially real). One must therefore complexify 
all fields, <p — > <f) R + irf* 1 . The Langevin equation then becomes 

a^R xc 

K R + V , K R =-Re^- , (2.3) 



6— ►<A K +i</ 



(2.4) 



Here we restricted ourselves to real noise. 

The complexification changes the dynamics substantially. Suppose that before 
complexification 6 is a variable with a compact domain, e.g. — 7r < 6 < it. After 
complexification, the domain is noncompact since — oo < 6 l < oo. Moreover, 
there will be unstable directions along which 6 l — > ±oo. This is best seen in 
classical flow diagrams, in which the drift terms are plotted as a function of 
the degrees of freedom 6 R , 6 1 . In Fig. [T] we show an example of a classical flow 
diagram in the XY model at finite chemical potential (to be discussed below). 
More examples of classical flow diagrams with unstable directions can be found 
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in Refs. [SI [TTJ [131 HSj- The arrows denote the drift terms (K R , K ) at (0 R , 1 ). 
The length of the arrows is normalized for clarity. In this case there are unstable 
directions at R ~ —0.7 and </> R ~ 2.4. The black dots denote classical fixed 
points where the drift terms vanish. Generally speaking, the forces are larger 
when one is further away from the fixed points. In absence of the noise, one finds 
generically that configurations reach infinity in a finite time, since the forces grow 
exponentially for large imaginary field values. 

When a Langevin trajectory makes a large excursion into imaginary direc- 
tions, for instance by coming close to an unstable direction, sufficient care in the 
numerical integration of the Langevin equations is mandatory. In some cases it 
suffices to employ a small stepsize e, after discretizing Langevin time as $ = ne. 
However, this does not solve instabilities in all situations. Moreover, a small 
stepsize will result in a slow evolution, requiring many updates to explore config- 
uration space. 

To cure both problems, we introduce an adaptive stepsize, e„, in the dis- 
cretized Langevin equations, 

R (n + 1) = <^{n) + e n Kf(n) + ^ Vx (n), (2.5) 
<P l x {n + l)=<Pl{n) + e n Kl{n), (2.6) 

where the noise satisfies 

(Vx{n)) = 0, (Vx(n)r) x '(n')) = 25 xx >5 nn t. (2.7) 

The magnitude of the stepsize is determined by controlling the distance a single 
update makes in the configuration space. Here we present two specific algorithms 
to do this. 

In the first formulation, we monitor, at each discrete Langevin time n, the 
quantity 

K™ x = max \K x {n) \ = max J Kf{n) + K x 2 {n). (2.8) 

x x V 

We then place an upper bound on the product eK max and define the stepsize e n 
as 

_(K max ) 



I/max 

n 



(2.9) 



Here e is the desired average stepsize (which can be controlled) and the expecta- 
tion value of the maximum drift term (K max ) is either precomputed, or computed 
during the thermalisation phase (with an initial guess). In this way, the stepsize 
is completely local in Langevin time and becomes smaller when the drift term is 
large (e.g. close to an instability) and larger when it is safe to do so. 

In the second formulation, we keep eK m3iX within a factor p relative to a 
reference value JC, i.e., 

-K < eK max <pJC. (2.10) 
p 
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If this range is exceeded the stepsize is increased/reduced by the factor p. This is 
iterated several times, if necessary. Both p and JC have to be chosen beforehand, 
but this does not require fine tuning as long as clearly inadequate regions are 
avoided. 

In the next sections, we apply these formulations to the XY model at nonzero 
chemical potential and QCD in the heavy dense limit respectively. 



3 XY model 

We demonstrate the first implementation using the three-dimensional XY model 
at finite chemical potential [17]. The action is 

2 

S = -f3 ^2 ^2 cos {(j) x - <f) x+ - i/jS ufi ) . (3.1) 

x u=0 

The theory is defined on a lattice of size Q = N T N^, with periodic boundary 
conditions in all three directions. The chemical potential is introduced as an 
imaginary constant vector field in the temporal direction [IS] and couples to the 
conserved Noether charge associated with the global symmetry <p x — > <p x + a. As 
always, the action is complex when /i ^ and satisfies S*(fi) = 5(— 
The drift terms, after complexification, read 

K ? = -p Yl [ sin 0* - A*) cosh W - - 

V 

+ sin (<f% - 4> x _ ) cosh {<p\ - (j) l x _ + /i5„ : o) , (3.2) 

K l = -i 3 J2 [ cos ~ ^+») sinh W ~ &+* ~ 



cos 



^_,)sinh(^,-0U + ^,o) • (3.3) 



As anticipated, they are unbounded due to the <p l variables. 

To construct the flow diagram in Fig. we have chosen the field variables at 
the six sites neighbouring <p x as random variables between ±7r. Note that the drift 
terms change sign when <p x — > <p x + ir, for given x, explaining the symmetry in 
Fig. [TJ The normalized drift terms and the classical fixed points (K^ = K\ = 0) 
are independent of (3. 

In an attempt to solve these Langevin equations numerically with a fixed step- 
size, we found that instabilities and runaway trajectories appear so frequent, that 



1 At nonzero chemical potential, it is preferable to interpret this system as a three- 
dimensional euclidean quantum field theory at finite temperature (with coupling (3 and in- 
verse temperature N T ), rather than as a three-dimensional classical spin system with inverse 
temperature /3. Models in this class can also be studied using world-line formulations [T9l 117) . 
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it is practically impossible to construct a thermalized configuration, even when 
the stepsize is very small, say, e ~ 10~ 5 . This becomes worse on larger volumes|§ 
We therefore switch to the adaptive scheme, using the first implementation. In 
Fig. |2]we show examples of the maximal drift term K ma,x and the adaptive step- 
size as a function of Langevin time for three different lattice volumes and two 
values of (3 and fi. We observe that the maximal force fluctuates over several 
orders of magnitude during the evolution (note the vertical logarithmic scale). 
Moreover, the frequency and size of the fluctuations increase when increasing the 
lattice volume. This is consistent with the picture developed above: on a larger 
volume there are more opportunities to be on a potentially unstable trajectory 
and subject to large forces. The effect also gets worse at larger chemical potential. 
i^ max and e n are inversely proportional, as expected. We note that although it is 
necessary to use a tiny stepsize occasionally, the algorithm is designed such that 
the evolution will continue with a larger stepsize as soon as possible. As a result, 
the time average of e n is close to the input timestep e = 0.01 in all casesjf] We 
emphasize that after the implementation of this algorithm we have not observed 
any instabilities, for a wide range of parameters (0.1 < /3 < 2, < fi < 6), 
lattice sizes (up to 16 3 ), and long runtimes (we explored millions of timesteps, 
corresponding to Langevin times of several thousand). 

To illustrate the method, we introduce two related models with a real action: 
the XY model with imaginary chemical potential fi = ifii, with the action 

Simag = -P/, COS (<t>x ~ <Px+0 + (J>l8vfl) , (3.4) 

and the phase quenched theory, obtained by taking the absolute value of the 
complex weight, i.e. e~ s —> \e~ s \ = e _5pq , which yields the action 

Spq = -fi'Y] cos (<j> x - (j) x+0 ) cosh (ji6 vfi ) . (3.5) 

X , V 

This is the anisotropic XY model, with direction-dependent coupling /3 V , where 
Po = /3cosh /j, and /3, = (3 (i = 1, 2). Both models are solved using real Langevin 
dynamics. The drift terms are bounded and there are no instabilities. 

In Fig. [3] we show the expectation value of the action density (S)/Q in the 
high-/? phase, at j5 = 0.55, for small values of The lattice sizes are 6 3 and 
8 3 , showing that finite size effects are under control. The result at imaginary 
H appears at ji 2 < 0, while the full and phase quenched results are plotted at 

2 This is in sharp contrast with the relativistic Bose gas in Rcf. j7j, where instabilities were 
not encountered for the parameter values used there. 

3 In practice, one may take e n < e always, to prevent the appearance of large timesteps. 

4 Rccall that at zero chemical potential, the three-dimensional XY model has a continuous 
phase transition at p c (/i = 0) « 0.4542 (see e.g. Ref. [2Q]), separating the symmetric phase at 
small P from the symmetry broken phase at large j3. 
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p=0.55. n=0.2S, 4 
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P=0.1,|l=2, 4 
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P=0.55, n=0.25, 16 




Figure 2: Example of the Langevin evolution of the maximal drift term K max / (3 
and the adaptive stepsize e n in the three-dimensional XY model with j3 = 0.55 
and /i = 0.25 (left) and — 0.1 and fi = 2 (right) on lattices of size 4 3 (top), 
10 3 (middle) and 16 3 (bottom). The input stepsize is e = 0.01. Note the vertical 
logarithmic scale. 
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Figure 3: Action density (S)/Q in the three-dimensional XY model as a function 
of fi 2 at /3 = 0.55 on lattices of size 6 3 and 8 3 , for the full theory (circle, square, 
/i 2 > 0), at imaginary \i (circle, square, /i 2 < 0), and phase quenched (triangles, 
/i 2 > 0). The vertical lines at \x\ = ir/N T indicate the Roberge- Weiss lines at 
imaginary \i. The dashed lines are the second-order fits (13.61 13.71) . incorporating 
the RW reflection symmetry. 

/i 2 > 0. At fi = all results agree (within the statistical error). The full and phase 
quenched results differ, as can be expected from e.g. a Taylor expansion of the 
observable for small /x. The results for imaginary and real /x appear continuous 
around /i 2 = 0, which is expected from the analyticity of the partition function 
in /i 2 on a finite lattice. The lines indicate second-order fits to the data on the 
6 3 lattice, with the results 

(S)/Q = -0.9433(7) - 0.502(4)/x 2 + 0.19(l)/i 4 , (3.6) 
(S) m /n = -0.940(2) -0.35(2)/x 2 - 0.04(3)//. (3.7) 

In the first case, the data at real and imaginary \x are combined in the fit. 

For imaginary /x we observe a cusp at \i\ = ir/N T . This is similar to the 
Roberge- Weiss transition in QCD [211 [22] and is due to the periodicity ji\ — > 
Note that the dashed lines in Fig. [3] reflect this symmetry. It would 

5 One way to see this is by using a field redefinition, <\> x = (f>' x + fi\T, which moves the \i 
dependence to the boundary condition <j>' N x = 4>' Q x — l*iiN T , similar as in fermionic models. 
The center symmetry is of course trivial in this model. 
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therefore be interesting to determine the phase structure of the XY model at 
imaginary chemical potential and finite N T . In the three-dimensional thermo- 
dynamic limit (N T is taken to infinity as well) and vanishing chemical poten- 
tial, the magnetized and symmetric phase are separated at the critical coupling 
/3 c (/i = 0). Consequently it is intriguing to analyse the interplay between the 
putative Roberge- Weiss transition and the standard magnetization transition in 
this limit, in particular since it would result in a breakdown of analyticity of 
Pdl^ 2 ) around /i 2 = and make /3 c (/i = 0) a multicritical point. 

4 Heavy dense limit of QCD 

To show the generality of the adaptive stepsize method we now apply it to the 
heavy dense limit of QCD in four dimensions. Here we shall present results 
obtained with the second algorithm. The stochastic quantization for this theory 
was studied in Ref. [6]. Here we briefly repeat the essential equations; we refer 
to Ref. [6] for further details. 

The gluonic part of the action is the standard Wilson SU(3) action, 

s B [u] = -p J2 E (l [ Tr u *#» + Tr U *>U - x ) ' t 4 - 1 ) 

fJ,<U 

where U X)llv are the plaquettes and (3 = 6/g 2 . The fermion determinant (starting 
from Wilson fermions) is approximated as 

det M « Y[ det (l + he^ T V x ) 2 det (l + he'^T^ 1 ) 2 , (4.2) 

X 

where h = (2k) Nt . Here k is the Wilson hopping parameter and A^ r = 1/T 
the number of sites in the temporal direction. The lattice spacing a = 1. The 
determinant refers to colour space only. The (conjugate) Polyakov loops are 

N T -1 

^ = n u ™+ k 1 = n ^x),4- (4.3) 

T=0 T = N T -1 

A formal derivation of Eq. (14.21) follows by considering the heavy (k — > 0) and 
dense (fi — > oo) limit, keeping the product /te M fixed (see Ref. [23] and references 
therein). The anti-quark contribution is kept to preserve the symmetry under 
complex conjugation. 

The Langevin process is 

U' x ,ti = R x,n u x,fi, Rx,n = exp [i\ a (eK Xfia + Ver] Xfia )] , (4.4) 
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Langevin iteration Langevin iteration 



Figure 4: Example of the Langevin evolution of the maximal drift term eK max 
(left) and the adaptive stepsize e (right) in the heavy dense limit of QCD with 
— 5, k — 0.12 and \x = 0.7 on a lattice of size 2 4 , using /C = 2 x 10~ 4 . 
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Figure 5: As above for ^Tr U4UI, indicating the deviation from unitarity during 
the evolution. 

where A a are the Gell-Mann matrices and the noise satisfies (r] xiia ) = 0, (r] Xfia Vy^b) = 
2S/ivS a b5xy The drift term K x/ia = -D x/M (S B + S F ), where S F = -IndetM, is 
complex due to the fermion contribution. For explicit expressions, see Ref. [6]. 

This Langevin process suffers from instabilities, which can be partially con- 
trolled using small enough stepsize. Applying the second formulation of the 
adaptive stepsize algorithm completely eliminates runaways. In this setup mea- 
surements are performed after a Langevin time interval of length A$ (this defines 
one "iteration"). The number of sweeps N sweep in an iteration depends on the 
Langevin timestep e: if the latter is decreased by a factor p, N swecp is increased by 
the same factor, and conversely. This makes the statistical analysis straightfor- 
ward, since each iteration has the same weight A& = eN swecp . Note that N Bweep 
will no longer be decreased and correspondingly e will no longer be increased if 
the former reaches 1. 
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To demonstrate this approach, we show in Fig. H] a characteristic evolution 
of eK ma * (left) and the stepsize e (right), using K = 2 x 10~ 4 and p = 2, with 
initial e = 1 x 10~ 5 , N swccp = 16. The total number of iterations is iVi ter = 3 x 10 5 
leading to a total Langevin time $ t ot = iViter x A?9 = 48. Varying these input 
parameters by factors of two or more does not change the results but may only 
affect the statistics. K ma,x is the maximum value of l-fQ^al over x, \i and a in the 
last sweep of an iteration. We observe that the product eK ma,x remains bounded, 
as required, while the stepsize (and hence K mBX ) fluctuate substantially. In Fig. [5] 
we show the evolution of |Tr U^Ul, which measures the deviation from unitarity 
[B], acknowledging the typical fluctuations in stationary regime. 

As mentioned above, after the implementation of this algorithm, we have not 
encountered any instability. 

5 Conclusion 

Instabilities can make complex Langevin simulations extremely problematic. We 
have shown that they result from large fluctuations in the magnitude of the drift 
term and the presence of unstable directions in complexified field space. This 
can be cured by using an adaptive time-local stepsize. The scheme is generic and 
could be applied to other theories such as QCD. While we have no proof that 
stable evolution is guaranteed, we have not encountered any instability using this 
method in the XY model at finite chemical potential and QCD in the heavy 
dense limit, for a wide selection of parameter values and system sizes. With a 
fixed stepsize on the other hand, instabilities can be so frequent that it is virtually 
impossible to generate a thermalized ensemble. 

Runaways are due to specific instabilities of the complexified Langevin equa- 
tions caused by the strong increase of the drift in the non-compact imaginary 
directions. We have shown that these runaways can be eliminated by using a 
dynamical step size, which indicates that they are not a question of principle for 
complex Langevin dynamics but one of numerical accuracy in following the tra- 
jectories. We consider this to be an important result of our analysis. Of course, 
for practical calculations one should consider further optimization of the algo- 
rithms by applying methods developed for general real Langevin processes (see, 
e.g., Ref. [21]) after analysing their adequacy to the complex Langevin problems 
of interest. 

We emphasize that the adaptive stepsize permits a fine tracing of the drift 
trajectories. If the process picks up a diverging trajectory the noise term will typ- 
ically kick the process off it. The present results suggest therefore that runaways 
are not due to following diverging trajectories j2Sj but rather due to following 
"wrong" trajectories, i.e. trajectories which, because of accumulating errors in 
the evaluation of the drift, do not belong to the dynamics of the problem. This 
both stresses the necessity of ensuring precision in the calculation and helps in 
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disentangling various sources of deficiency in the search for a reliable method. 
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